Self-adapting method for the locahzation of quantum critical points using Quantum 

Monte Carlo techniques 
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A generalization to the quantum case of a recently introduced algorithm (Y. Tomita and Y. Okabe, 
Phys. Rev. Lett. 86, 572 (2001)) for the determination of the critical temperature of classical 
spin models is proposed. We describe a simple method to automatically locate critical points in 
(Quantum) Monte Carlo simulations. The algorithm assumes the existence of a finite correlation 
length in at least one of the two phases surrounding the quantum critical point. We illustrate 
these ideas on the example of the critical inter-chain coupling for which coupled antiferromagnetic 
S = 1 spin chains order at T = 0. Finite-size scaling relations are used to determine the exponents, 
V = 0.72(2) and r) — 0.038(3) in agreement with previous estimates. 
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For studies of phenomena related to phase transitions, 
either classical or quantum, a precise determination of 
the critical point is of crucial importance, especially so 
in numerical work. For simulations of classical statistical 
models several algorithms [1,2] have been proposed, that 
automatically adjust to the critical point. Using ideas 
from invasion percolation Machta et al. [1] proposed an 
invaded cluster algorithm where critical clusters are gen- 
erated directly, yielding the critical temperature as out- 
put. More recently, Tomita and Okabe [2] have suggested 
using a simpler method called probability-changing clus- 
ter algorithm to localize the transition temperature of 
Ising and Potts models within the framework of clus- 
ter algorithms. In Swendsen and Wang or Wolff algo- 
rithms [3], they propose to, during the simulation, either 
increase or decrease the probability p of connecting spins 
in the same cluster, and hence the temperature, accord- 
ing to whether the last clusters were percolating or not. 
Thus, systematically moving towards the cluster distri- 
bution at the critical point. 

The generalization of either algorithm to the simula- 
tion of quantum systems would clearly be highly desirable 
but is difficult for several different reasons. Firstly, the 
structure of the clusters are more complicated than for 
classical systems and their relation to percolation physics 
is less clear. During the generation of a cluster it is not 
possible to assign a probability according to which a sin- 
gle "spin" is added to the cluster. Secondly, the dynam- 
ical critical exponent z, that links real and imaginary 
time correlations, is usually not known [4]. Hence, a di- 
rect generation of critical clusters would appear difficult 
due to the lack of knowledge of the relation between real- 
space and imaginary time correlations. A resolution of 
these difficulties would be very interesting and here we 
discuss a first step in this direction. We propose to gen- 
eralize the approach of Tomita and Okabe [2] to Quan- 
tum Monte Carlo (QMC) simulations performed using 



the loop algorithm [5]. We show that adjusting the criti- 
cal parameter during the simulation according to whether 
an estimate of the correlation length is larger or smaller 
than some threshold, allows for a very precise determi- 
nation of the quantum critical point. With relatively 
modest effort most of the associated critical exponents 
can also be determined. 

The method we propose is closely related to the one 
suggested by Tomita and Okabe [2]. During the simula- 
tion a variable measuring the criticality of the system, in 
our case the correlation length ^, is continuously moni- 
tored and one then changes a control parameter, A, ac- 
cording whether or not this critical variable satisfies a 
criterion. Here A can be the temperature fcsT = as 
in Ref. [2] or directly a parameter of the Hamiltonian, 
controlling the transition. The algorithm should then 
equilibrate A to an average critical value Ac(L), depend- 
ing on the system size L, where the specified criterion is 
satisfied as an equality. The most difficult thing one has 
to do with this algorithm is to select a simple numerical 
criterion indicating whether the critical variable is in the 
ordered phase or not. As discussed in [2], by determin- 
ing the critical value Xc{L) for different linear sizes L of 
systems, it is possible to obtain Ac by extrapolating to 
the thermodynamic limit. 

Let us briefly describe the method for a system of lin- 
ear size L. First, one has to choose an initial value Ao(i) 
for the parameter A. Ao(i) can be chosen according to 
a physical intuition of the critical point, or as the result 
of a very short preceding simulation at this size, or of 
the simulation for a shorter linear size. In most cases 
it is useful to choose Ao(i) to be in a phase with a fl- 
nite correlation length. Then a certain number iVwarm of 
"warmup" Monte Carlo steps are made to help the sys- 
tem to locate approximately Ac(i). During these steps, A 
is modified by a certain amount ±AoA every Nt sweeps, 
according to our criterion which we will discuss later. 
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The Nt sweeps are needed to "thermalize" the system 
to the current value of A to avoid autocorrelation prob- 
lems (if A is changed at every step, the configurations 
might be too correlated). Note that AqA is not necessar- 
ily fixed, and can for example be changed linearly from 
AqA to the smaller fixed value AA used in the ensuing 
part of the calculation as discussed in [2]: this permits 
the system to take big steps in A in the early stages of the 
simulation if Xo{L) is far from the critical value Xc{L). 
When the A'warm steps are finished, one makes A^mcas 
steps of measurement with the same procedure : every 
Nt steps, A is changed by a small and fixed value ±AA 
and estimates of observables (including A) during the Nt 
steps are recorded. At the end of these iVmeas sweeps, we 
get a distribution of A centered around the critical value, 
Xc{L), for this system size. For the interpretation of the 
final distribution of A it is important to have used a fixed 
value of AA during the A^meas steps of measurement. We 
usually choose AA ^ where V is the volume of the sys- 
tem. Finally, in order to determine Ac, the calculation is 
repeated for different linear sizes L and an extrapolation 
to the thermodynamic limit is performed, in principle 
allowing for a determination of the correlation length ex- 
ponent v. Other critical exponents can be determined in 
this step with the scaling of other observables. 

Now we turn to a discussion of the criteria used to lower 
or raise A. In Refs. [1,2] two criteria based on the perco- 
lation of the clusters were used. These criteria, which are 
referred to as the extension and topological rule in the 
related invasion cluster algorithm proposed by Machta et 
al. [1] for classical spin systems, can in principle also be 
used in QMC simulations with the loop algorithm: per- 
colation or winding of loops can be used to determine 
the change in A. However, the interpretation is in this 
case less clear. Most notably the exponential distribu- 
tion of the size of the clusters is very wide and not very 
useful for a determination of Ac. In particular, due to 
finite-size effects, it is not likely to change qualitatively 
when A is tuned through the transition. Eventually, an 
algorithm based on the cluster distribution in the ther- 
modynamic limit at T = would be interesting to pur- 
sue, but the present method is in our opinion simpler and 
more general. Hence, we choose a different criterion: A 
is modified according to an estimate of the correlation 
length ^ during the Nt steps. It is the same kind of crite- 
rion used in the fixed cluster algorithm [6]. It should be 
noted that this criterion is not necessarily related to the 
loops of the algorithm and can probably be used in other 
Monte Carlo simulations which do not admit a cluster al- 
gorithm. However, here the loop algorithm is very helpful 
because it quickly equilibrates configurations during the 
Nt sweeps and because it allows the use of improved es- 
timators [7] for correlation lengths. Note the importance 
of using a reasonably large Nt to have a good estimate 
of the correlation length. 

In the following, we will apply and check this method 



to quantum S = 1 spin systems in order to determine 
the critical value of the perpendicular coupling, J-^ , nec- 
essary for antiferromagnetic order at T=0 for coupled 
S = I spin chains. The S = 1 spin chain is known to be 
disordered [8] with a gap Ah — 0.41J and a finite cor- 
relation length ^ ~ 6 [9]. On the other hand, the S = 1 
square antiferromagnet is Nccl ordered at T = [10] 
and it is known that a finite inter-chain coupling J^, at 
which the system orders, exists [11-13]. We consider the 
following Hamiltonian: 

with integer S = 1 spins and periodic boundary condi- 
tions. Previous work on this transition [11-13] has esti- 
mated the critical coupling and the exponents with the 
most recent estimate [13] for being = 0.043648(8), 
a value surprisingly small compared to the Haldane gap 
Ah ^ 0.41J. 

We investigate this quantum phase transition by set- 
ting our parameter A as being the inter-chain coupling 
J-"-. Our first criterion (A) is the following: if the cal- 
culated correlation length along the chains is larger than 
the system size (^a, > L), then the system is in an or- 
dered phase for this value of and we decrease J-"-. If 

< L,we increase J-"-. The dividing criterion is thus: 

A: C = i. (2) 

We then record the distribution of J^, and estimate 
J^{L) for each lattice size. This criterion gives a good 
estimate of the value of the critical point, but will fail in 
determining it exactly and in determining critical expo- 
nents. This is due to the fact that an estimate of ^ for a 
system size where ^ ~ L suffers from pronounced finite- 
size effects. In particular, the extrapolation of J^iL) to 
the thermodynamic limit does not follow a simple scaling 
law allowing for a determination of the critical exponents. 
It is a well-known numerical fact that finite size correc- 
tions are significantly reduced when the lattice size is a 
few times larger than the calculated correlation length. 
Typically, if 6^ < L, finite-size corrections should be rel- 
atively unimportant. Hence, we use a second criterion 
(B) 



in order to reduce the finite-size corrections. 

We use the continuous time QMC loop algorithm [5] 
where spins 1 are simulated by two symmetrized spins 
1/2 [14] in the imaginary time direction. In order to 
study this quantum phase transition it is crucial to set 
T ~ and for this purpose we mainly use the recently 
proposed improvement of loop algorithm [15] which al- 
lows for simulations directly at /3 = oo. Observables are 
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measured with the help of improved estimators [7]. In 
particular, the correlation length ^ is measured with the 
second moment method [16,7,14]. In general, simulations 
have been made on LxLj, lattices. For the smaller lattice 
sizes, square lattices were used with T = [15]; for the 
largest lattice sizes Ly < L (with Ly » ^y) was used 
and simulations were performed at a finite temperature 
significantly smaller than the smallest gap in the system. 
We typically set Nt = 10^, N^^,^ = 10^ and TV^eas = 10^ 
sweeps. 

We now turn to a discussion of our results. In Fig. 1 
we show an example of the distribution of J± obtained 
during the simulation for L = 66 using criterion B. The 



0.04 



0.03 



0.02 



0.01 



0.025 



0.05 



0.075 



0.1 



0.0232 




0.024 



figure 1: Histogram of obtained for the L = 66, Ly = 
10 lattice at T = with the criterion B = The 
distribution is peaked around J^(66) = 0.0236(1). 



distribution is sharply peaked around a well defined mean 
value allowing for an easy determination of J;^(L). In 
principle additional information pertaining to the critical 
exponents should be extractable from a detailed study of 
the scaling of the complete distribution with the system 
size L. For the simulations we have performed we did 
not have sufficient statistics to exploit this possibility. In 
Fig. 2 we show J^{L) as a function of L for the two cri- 
teria, A and B. As previously mentioned the criterion A 
suffers from pronounced finite-size effects. In particular, 
criterion A always over estimates the critical coupling. 
This might appear counter intuitive since the correlation 
length is expected to diverge at J^. However, due to 
finite-size effects, the estimated correlation length never 
surpasses the linear size of the system until well inside 
the antiferromagnetic phase. One should note that for 
a finite quantum system there is always a finite gap re- 
straining the correlation length. In order to extract the 
critical exponent we can as a first approximation fit the 
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figure 2: J^{L) as a function of 1/L for the two different 
criteria A (top), B (bottom) used. Lines are fits to the 
equation (4) and give = 0.0428(4) for the criterion 
A and = 0.0438(4) for the B criterion. The values 
of Jc{L) found for the first (second) criterion are always 
larger (smaller) than the critical coupling. 



results obtained to a power law: 

Jt{L) = Jt+AL-\ 



(4) 



Here is the critical coupling (in the thermodynamic 
limit) and A a constant. For the criterion B (^^ = L/6), 
a = l/v where v is the correlation length exponent. For 
the criterion A, a ^ 1/v due to finite-size corrections 
and the L dependence is non-trivial. We obtain = 
0.0438(4) and v = 0.77(2) for the criterion B, and = 
0.0428(4) for the criterion A. The more reliable values are 
of course the first ones, but one can note that even with 
the criterion A, we have a good estimate of J^. We can 
try to improve on this estimate by including corrections 
to scaling. In analogy with previous work on finite-size 
scaling [17] we include a first order finite-size correction 
to the correlation length. 



e = a{J^ - J^{L)y 



(5) 



Including this finite size correction we find v = 0.72(2), 
= 0.0436(2). In essence this correction amounts to 
a relatively precise estimate of the correlation length ^ 
in the thermodynamic limit for the given value of J^{L). 
This is then fitted to the expected form for the divergence 
of the correlation length as is approached. In princi- 
ple higher order corrections could be included, we have 
verified that in the present case they are extremely small. 
This value of is in agreement with the recent QMC 
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"traditional" simulations [13] and is more precise than 
older results obtained with exact diagonalizations [11] or 
QMC [12] calculations. 

Once a reliable estimate of J^iL) is obtained high- 
precision estimates of additional critical exponents can 
be obtained by performing calculations directly at J^{L). 
In Fig. 3 we show results 
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save a lot of computational effort. In particular, no a pri- 
ori knowledge of z is needed. Most other methods, as for 
instance the use of Binder cumulants [21], would require 
the estimation of a critical quantity for many different 
system sizes for each value of the control parameter. The 
present method only requires the knowledge of the crit- 
ical quantity for one system size for each of the control 
parameter. As mentioned in the introduction, it would 
be highly desirable to extend also the invaded cluster al- 
gorithm [1] to the simulation of quantum systems and 
we hope the present work could be a first step in that 
direction. 
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figure 3: The staggered susceptibility Xs as a function of 
L on a log-log scale. The power-law fit (6) displayed as 
the line yields r? = 0.038(3). 

of the size dependence (on log-log scale) of the stag- 
gered susceptibility Xs with the criterion B (^x = L/Q). 
According to standard finite-size scaling theory this 
quantity should scale as 
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(6) 



even with J^{L) a function of L. Additional finite size 
corrections of the form (5) for Xs were not included. Fit- 
ting to this form (6), we obtain r] = 0.038(3) in complete 
agreement with previous estimates. The obtained values 
of rj and v for the critical exponents are in agreement 
with Matsumoto et al.'s results [13]. The exponents are 
very close to the values for the 3D classical Heisenberg 
model [18] confirming the expectation that two dimen- 
sional quantum phases transitions of quantum Heisen- 
berg magnets belong to the university class of classical 
magnets [19] as seen in in other simulations [20,13]. 

In this paper, we have proposed a very simple method 
to determine critical points in QMC simulations as a 
generalization of the probability-changing cluster algo- 
rithm [2] and applied it to determine the critical value 
of inter- chain coupling for spin S = 1 chains and calcu- 
late critical exponents for this quantum phase transition. 
The results are in agreement with previous estimates. It 
is important to underline that this method is easy to em- 
ploy and adapt to other phase transitions and that it can 
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